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Abstract. Boolean threshold networks have recently been proposed as 
useful tools to model the dynamics of genetic regulatory networks, and 
have been successfully applied to describe the cell cycles of S. cerevisiae 
and S. pombe. Threshold networks assume that gene regulation processes 
are additive. This, however, contrasts with the mechanism proposed by 
S. Kauffman in which each of the logic functions must be carefully con- 
structed to accurately take into account the combinatorial nature of gene 
regulation. While Kauffman Boolean networks have been extensively 
studied and proved to have the necessary properties required for model- 
ing the fundamental characteristics of genetic regulatory networks, not 
much is known about the essential properties of threshold networks. Here 
we study the dynamical properties of these networks with different con- 
nectivities, activator-repressor proportions, activator-repressor strengths 
and different thresholds. Special attention is paid to the way in which 
the threshold value affects the dynamical regime in which the network 
operates and the structure of the attractor landscape. We find that only 
for a very restricted set of parameters, these networks show dynamical 
properties consistent with what is observed in biological systems. The 
virtues of these properties and the possible problems related with the 
restrictions are discussed and related to earlier work that uses these kind 
of models. 

1 Introduction 

The analysis of the dynamics of genetic regulatory networks in living organisms 
is a complicated task and a central challenge in current research for a complete 
understanding of complex biological systems. Historically, the dynamical be- 
haviour of the biochemical elements in small genetic circuits has been accurately 
described using differential equations, which capture the underlying reaction- 
diffusion kinetics that take place in these systems |H2l3j . However, this approach 
faces important difficulties for the modeling of large genetic networks, being the 
main difficulty that these mathematical models may involve a very large amount 
of parameters. This is a serious problem both practically and theoretically. Prac- 
tically since these parameters may be largely unknown for many systems, and 
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theoretically since such a detailed description may obscure the essential proper- 
ties of the regulatory processes in the systems under consideration |4|5l . Because 
of this, Boolean networks have recently been increasingly used as the best hrst 
approach for the modeling and understanding of the essential properties of real 
regulatory systems that incorporate large amounts of data |6I7| . 

Boolean networks have been extensively studied for decades [5] , and were in- 
troduced for the modeling of large regulatory systems by S. Kauffman as a first 
attempt to understand the general dynamical properties of the gene regulation 
and cell differentiation processes |9ll0j . However, it was only recently that the 
necessary information to test them on real biological genetic networks has been 
available. Examples are models of the genetic network of flower development 
in Arabidopsis thaliana |llll2j . the regulatory network determining embryonic 
segmentation in Drosophila melanogaster |13| , the network controlling the differ- 
entiation process in Th cells [14] , the cell cycle networks of Saccharomyces cere- 
visiae |15j and Saccharomyces pombe |16j . among others. One of the advantages 
of the Boolean approach is that it is not necessary to know the kinetic details 
of the interactions (e.g. promoter affinities, degradation constants, translation 
rates, etc.). Rather, only the logic of the regulatory interactions is needed, such 
as the specific activatory or inhibitory nature of the genetic regulations [415] . By 
incorporating this information, available nowadays from high-throughput exper- 
iments, into the Boolean approach, it has been possible to predict the temporal 
sequence of gene activities as well as the stable and periodic patterns of gene 
expression in wild type and in many mutants of the organisms mentioned before. 

There are, however, important differences in the way in which Boolean mod- 
els have been implemented by different groups, and it is not clear whether or not 
these different implementations would yield equivalent results. The general for- 
mulation of the Boolean network model is the following. We assume that the net- 
work is represented by a set of N Boolean variables (or genes) {<7x, ^2, ■ ■ ■ , cat}, 
each of which can be in two different states cr, = 1 (active) and cr, = (inac- 
tive). The state of each gene <7j is controlled by fcj other genes of the network, 
{oij, . . . , <Ji k , }, which we will refer to as the regulators or the inputs of o~i. The 
number A:, of regulators of each gene depends on the topology of the network in 
such a way that the probability for a randomly selected node to have k regula- 
tors is given by the probability distribution Pj„(fe). Once every gene has been 
provided with a set of regulators, the dynamics of the network are given by the 
simultaneous updating of all the gene states according to 

o-i(t + 1) = F t ( ffll ((),(7, 2 (t), . . .,a ik . (i)) , (1) 

where Fi is a regulatory function, specific to the gene a^, that is constructed 
according to the activatory and inhibitory nature of the regulators of <7j. 

The differences in the implementation of the Boolean approach mentioned be- 
fore are related to the way in which the regulatory functions Fi are constructed^ 



1 There is another important difference in the Boolean implementation which is not 
related to the regulatory functions but that is worth mentioning, which is the syn- 
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For instance, the regulatory functions used in |ll|12j for the A. thaliana flower 
development network were carefully constructed taking into account the current 
biological knowledge about the combinatorial action of the regulators on their 
target genes. This combinatorial action takes place, for instance, with dual regu- 
lators whose inhibitory or activatory nature on the target gene depends upon the 
presence or absence of other regulators (which may compete for the same binding 
site in the promoter region) |17j . In contrasts, the regulatory functions used in 
[15] and |16j for the cell cycle networks of S. cerevisiae and S. pombe are thresh- 
old functions similar to the ones used in artifical neural networks [18119] . These 
two schemes, combinatorial functions vs. threshold functions, are very different 
not only mathematically, but in their very nature. For the use of threshold func- 
tions requires the strong assumption that the effect of activatory and inhibitory 
regulations, rather than combinatorial, is simply additive. 

In spite of this strong assumption, Boolean models with threshold functions 
seem to predict the correct biological sequence of events in the cell cycles of 
S. cerevisiae and S. pombe |15|16j . The dynamics in each of these systems ex- 
hibit one big attractor that corresponds to the experimentally observed stable 
state at the end of the cell cycle. This result suggests that, under certain condi- 
tions, gene regulatory interactions can indeed be considered as purely additive. 
In such cases, Boolean models with threshold functions are useful to describe 
real genetic networks and understand their dynamical properties [20] . There- 
fore, a thorough study of these kind of mathematical models is necessary. How- 
ever, although Boolean networks with threshold functions have been extensively 
studied in the context of spin glasses [H]-[2E] and artificial neural networks 
[18 19 , their dynamical properties in the context of gene regulation are largely 
unknown. For only the most simple cases of fixed connectivities, equal activa- 
tor/repressor sterngths and proportions, and fixed threshold values have been 
explored [H27]Q 

In this work we investigate the generic dynamical properties of Boolean net- 
works with threshold functions. Our main goal is to compare the behavior of 
these threshold networks with the one that is already known for standard ran- 
dom Boolean networks (also termed Kauffman networks), focusing on the prop- 
erties that are relevant to gene regulation processes. To this end, we use different 
connectivities, activator/repressor strengths and proportions, and threshold val- 
ues. In Sec. [2] we describe the Boolean threshold network model and present 
examples of strong deviations from the "normal" behavior observed in Kauff- 
man Boolean networks. Next, in Sec. [3] we use the annealed approximation [28] 
and the average influence [29] of these networks to calculate the phase diagram 



chronous versus the asynchronous updating schemes. Throughout this work we will 
use synchronous updating because we want to focus on the differences regarding the 
construction of the regulatory functions. 
2 Usually, for spin glasses and neural networks the nodes Oi take the values {+1, — 1} 
(rather than {1,0}). Although models using the spin-like values {+1,-1} can be 
mapped onto models using {1,0}, the mapping requires a fine tuning of the threshold 
values 8i. 
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for the different parameters involved. In Sec. [4] we present numerical evidence 
to support the analytical results and discuss the case where anomalies between 
the theoretical prediction and the numerical simulations arise. Finally, we dis- 
cuss and summarize our results, highlighting their implications in terms of the 
applicability of threshold networks for the modeling of gene regulation. 



2 The Boolean threshold network model 
2.1 Definition and general properties 

In what follows we will refer to Boolean networks with threshold functions as 
Boolean threshold networks (or BTN's). Since threshold functions are a subset 
of the general class B comprising all possible Boolean functions, it is clear that 
BTN's are a subset of the ensemble of random Boolean networks (RBN's) intro- 
duced by Kauffman, in which the regulatory functions -Fj are randomly chosen 
from B. In the context of gene regulation, the dynamics of BTN's are given by 



<n(t + l)=Fi (a ix (t), a^t), a ik . (i)) 
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where 



(2) 

| are the fcj regulators of <jj. The interaction strength (or 

weight) dij takes a positive (or negative) value if <7j is an activator (or a repres- 
sor) of <7j, respectively]^] The activation threshold 0$ of cr, indicates the minimum 
value of the sum required for the activation of the node to take place. The dy- 
namic rule given in Eq. |2]) is the same as the one used in Refs. [15116] . However, 
in that work the authors considered the simple case in which a^j = 1 for acti- 
vators, aij = —1 for repressors, and 6i — for almost all the nodes except by a 
few ones. Additionally, "self-degradation" was introduced to some of the nodes 
just by making a,,, = — 1. 

The number ki of regulators for each node o~i is drawn from a probability 

distribution Pj„(fc), and then these regulators iu^, ■ ■ ■ ,o~i h . j are randomly cho- 
sen from anywhere in the system. Each regulatory interaction strength a^j is 
set activatory with probability p and inhibitory with probability 1 — p. All ac- 
tivatory interactions have a value a^j = ag, whereas the inhibitory interactions 
have a value a^j = — or, where ac and an are positive integers. The ratio ac/cm 
measures the relative importance of activation over repression. Thus, if ac/aa 
is small, then inhibitory interactions are dominant, whereas if clc/cir is large, 



3 Ofcourse, ctij — if there is no interaction between ai j and <Ti. 
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then activation dominates over repression. Finally, we set a fixed value of the 
activation threshold 9i = 9 for all nodes. We consider three cases corresponding 
to three different threshold values: 9 = 0.5, 9 = and 9 = —0.5. The rationale 
for this choice is two fold. First, these values suffice to illustrate the effects of 
integer and non-integer thresholds on the dynamics. And second, because these 
are the values that have been used in models of real genetic networks, obtaining 
good agreement with experimental observations |15|16j . 

On thing that should be noted from Eq. ^ is the effect that the value of 
the threshold 9i has on the dynamics. If we consider only integer values for the 
interaction strenghts djj, then the equality in Eq. ^ can be attained only if $i 
is also an integer. In such a case, the last row on the right-hand side of Eq. ([2| 
implies that every node regulates itself. In other words, given the interaction 
strengths djj and the thresholds 9 i: the right-hand side of Eq. ^ can be writ- 
ten as a Boolean function Fi only if we assume that <jj belongs to its own set 
of regulators. Because of this, we have explicitely written <ii{i) as one of the ar- 
guments of the regulatory function Fi. This self-regulation does not necessarily 
happen in Kauffman Boolean networks, and it can make a big difference with 
regard to the dynamical behavior. As we will see below, the fact that integer 
threshold values allow the node <Xj to simply stay in their previous state and 
essentially freeze plays a mayor part in the dynamical behaviour of the network 
and in its use for biological modeling. 

Note also from Eq. ^ that all the information necessary for the network 
dynamics is contained in a A-dimensional vector = {6\, 62, • ■ ■ , On) whose 
components are the thresholds, and a N x N matrix A. This matrix is such 
that [A]ij = dij if <7j. is a regulator of o^, and [A]ij = otherwise. This is 
very different from what happens in RBN's, where to store all the information 
necessary for the network dynamics we need a N x N matrix containing the 
topology of the network, and a Boolean function Fi for each node crj. Each of 
these functions has 2 ki entrances, one for each configuration of its ki inputs. 
As mentioned before, for a given set of thresholds and interaction strengths, we 
can also create a Boolean function corresponding to the rule given in Eq. (§. 
However, this limits the set of possible Boolean functions that can be obtained. 



2.2 The Derrida map : Deviations from the Kauffman behaviour 

One of the most useful ways to study the general dynamics of Boolean net- 
works has been in terms of the propagation of perturbations (also called damage 
spreading) throughout the network. To this end, let us denote as E t the dynami- 
cal configuration of the network at time t, that is, Et = {cr\ (t) , 01 (t) , . . . , a jv(i)}. 
Let Eq and Eq be two slightly different intitial configurations, namely, Eq is al- 
most identical to E except by a few nodes which have reversed their values (this 
is the initial perturbation or the initial damage). Under the dynamics given in 
Eq. ([2]), each of these initial configurations will generate a trajectory throughout 
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time: 

Sq — > Si — » • • • — > Z't — » • • • 

£ -> £ ! -»■ >S t ^--- 

These two trajectories may eventually converge (the initial perturbation dis- 
appears), diverge (the initial perturbation amplifies), or remain "parallel" (the 
initial perturbation neither grows nor disappears). These three different behav- 
iors determine the dynamical regime in which the network operates: In the or- 
dered regime the two trajectories typically converge after a transient time. In the 
chaotic regime the system becomes very sensitive to small changes in the initial 
condition and the two trajectories diverge from each other. The intermediate 
case where, on average, perturbations retain their same size corresponds to the 
so called critical regime. The critical regime has been extensively studied and 
appears to be characteristic property of genetic networks [30 31 32 33|34| . 

We quantify the propagation of perturbations in the network in terms of the 
time evolution of normalized Hamming distance h(i), which is defined as 

, N 

h(t) = d(£ t ,Z t )=-J2\a i (t)-a i (t)\- (3) 

i=l 

The assymptotic value hoo = hm h(t) is the final size of the avalanche of per- 

t— >oo 

turbations and acts as the order parameter of the system: In the ordered regime 
hoo = 0, while in the chaotic regime hoo > 0. In the critical regime lim h(t) = 

t— >oo 

only marginally, which means that it can take a long time for a small perturba- 
tion to disappear. 

For a given network realization, hoo can be computed numerically in two 
different ways. The first way is a direct implementation of the definition. We 
start out the dynamics from two different initial conditions Sq and Eq, and let 
the system evolve for a long time t r . Then, /i M is the Hamming distance h(t r ) 
between the two final configurations £ tr and S trl averaged over many pairs of 
initial conditions. We will denote the value of the order parameter obtained by 
this method as h^. 

The second way to compute hoo is by means of the so-called Derrida map 
M (h) [35], which relates the size of a perturbation avalanche between two con- 
secutive time steps, that is, h(t + 1) = M (h(t)). Starting from two different 
initial configurations whose Hamming distance is ho, succesive iterations of this 
map eventually converge to hoo- Thus, hoo is the stable fixed point of the Der- 
rida map: hoo — M (hoo)- For RBN's, mean-field theory computations show that 
M (h) is a continuous convex monotically increasing function with the properties 
M(0) = and M(l) < 1. For threshold networks this mapping is still continuous 
and satisfies M(0) = and M(l) < 1, but it is not clear whether or not it is a 
monotonically increasing function. Nonetheless, for the set of parameters we use 
in this work M (h) seems to satisfy all the properties predicted by the mean-field 
theory. The fulfillment of these properties is important because this guarantees 
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the existence of one and only one stable fixed point. In this case, the dynamical 
regime in which the network operates is determined by the slope at the origin 
of M (h), called the average network sensitivity S: 



s = dM{h) 



dh 



(4) 

h=0 



If S < 1 then hoo = and the system is in the ordered phase, whereas if S > 1 
then ft-oo > and the system is in the chaotic regime. The critical regime is 
attained for S — 1, which is the point at which the phase transition between the 
ordered and chaotic regimes occur |35l36j . 

To compute M (h) numerically for a given network realization, we start from 
two different configurations S and S separated by a Hamming distance ho- 
Next, we evolve these two initial configurations just one time step and compute 
the Hamming distance hi between the resulting configurations Si and Si- The 
value M [ho] of the Derrida map at ho is then obtained by averaging hi over 
many pairs of initial conditions whose Hamming distance is ho- By doing this 
for all values of h £ (0, 1) we can construct the full curve M(h) and compute 
its fixed point h^ . We will denote the value of the order parameter obtained by 

(2) 

this method as hid- 

For general RBN's it has been shown that hoo and ftS are very close to 
each other. Actually, in the thermodynamic limit N — > oo they are the same 
{35] . The reason for this is that in RBN's the temporal correlations between two 
consecutive configurations S t and S t +i are inversely proportional to the num- 
ber of nodes N . Therefore, for large networks with completely random Boolean 
functions the mean-field conditions are satisfied and the temporal evolution is 
essentially dependent on the previous time step only. However, when temporal 
correlations extend over several time steps, the Derrida map does not accurately 
predict the value of the order parameter. In such cases h^ and h& can differ by 
a large ammount. This non-ergodic behavior in the network dynamics has been 
observed in Boolean networks in which only a small subset of the class of all 
Boolean functions are used [37] , such as canalyzing functions 29J and threshold 
functions with equal values and proportions of activation/repression strengths 
[27] . For the general case of RTN's we also observe a large deviation from the 
ergodic behavior assumed by the mean-field computation. 

In Fig.[l]we plot hiV (diamonds) and h^} (circles) as functions of the network 
connectivity K, for RBN's (Fig. [lji) and RTN's (Fig. \fy,c). We also plot the 
quantity h*^ predicted analytically using the annealed approximation presented 
m Sec.§ which is a generalization of the one reported in Ref . [27] Q Note from 
Fig. [IJi that for RBN's the three values of hoo are identical within numerical 
accuracy, which reflects the ergodicity of the system in this case. However, for 



RTN's such ergodicity dissappears, as it is apparent from the fact that hoo is 
quite different from hW ■ Fig. lib corresponds to the case in which 9 = for all 



4 This computation incorporates in an approximate way the temporal correlations 
between succesive network states using the final number of active and inactive nodes. 
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Fig. 1: The nonergodicity of the system is illustrated by plotting the order pa- 
rameter h& computed directly from the definition (red diamonds) , and the order 

(2) 

parameter haJ computed as the fixed point of the Derrida map (blue circles). 
The analytic prediction h*^ from the annealed computation presented in see 
Sec.|3]is also plotted (solid line). Three different ensembles of networks are used: 
(a) Standard Kauffman networks (RBN's). In this case — /loo = h^, which 
shows that RBN's are ergodic. (b) Random Thrcsdhold Networks (RTN's) with 
p = 0.5, clq = an , = 1 and 9 = 0. Note in this case that hx ^ hx (although 

(2) 

hcJ = which reflects the nonergodicity of the dynamics. Finally, (c) cor- 

responds to RTN's with p = 0.5, ac = 1, clr = 20 and 9 = 0. In this last case, 
not only is hx ^ , but also the analytic prediction h*^ fails completely. In 
all cases, each point is the average over 100 network realizations, each having 
N = 1000 nodes. For each of these networks we used 10000 pairs of random 
initial conditions. 



nodes and the weights take the values ajj = ±1, chosen with equal probability. 
This strong deviation is surprising, especially since it happens for the simplest 
case similarly to the one used in Refs. |15|16j for the modelling of the yeast 
cell-cycle networks. Furthermore, departure from ergodicity is even worse for 
unequal weights, as shown in Fig. [T];, where the negative weights were chosen 
to be ten times stronger than the positive weights, i.e. = 10 and clq = 1. In 
this case h& does not only deviates from , but also the analytical prediction 
h* x completely fails to reproduce /loo ■ The above results indicate that some care 
must be taken when choosing the parameters in RTN's if these networks are 
to be used for biological modelling. We explore this issue furtherly in the next 
sections. 
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3 The phase diagram 



The Derrida map M (h) can be computed analytically within the context of the 
so-called annealed approximation, first introduced by Derrida and Pomeau |28j . 
This mean-field technique assumes statistical independence between the nodes 
and neglects the temporal correlations developed throughout time between suc- 
cesive states of the network. The annealed approximation has been successfully 
used in RBN's to obtain analytically where the phase transition occurs for differ- 
However, the mean-field assump- 



ent topologies and network parameters |36|8j 
tions fail dramatically for RTN's as it is illustrated in Fig. [T] In an attempt 
to improve the annealed approximation, one has to incorporate into the analy- 
sis the temporal correlations between succesive network states |29l37j . This has 
been done for particular values of the parameters |27) . Here we present a gen- 
eralization of this computation valid for different activator/repressor strengths, 
ratios and thresholds. 

We start the computation of M(h) by introducing the quantity I^ kd \ known 
as the influence of kd variables. Let us consider an arbitrary network in the 
ensemble of RTN's, and pick out a node <7j with fc, inputs. Let S t and S t 
be two network configurations in which kd of the inputs of o~i (with kd < k{) 
have been damaged, namely, these kd inputs have opposite values in these two 
configurations^] /( fc<J ) is defined as the probability that this initial damage of kd 
inputs propagates one time step, which means that o~i will have different values 
in St+i and ^t+i- These influences do not only depend on the ensemble of 
Boolean functions used, but have also been shown to depend heavily on the bias 
in the expected probability with which the system visits the different states of 
its configuration space. In previous work, this bias has been expressed in terms 
of the fraction b(t) of active nodes in the system |27|38)37| . 



By using the annealed approximation assumptions, in Appendix A we show 
that the temporal evolution of b(t) is given by 



b(t + 1)=B (b(t)) = p+ (b(t)) + b(t) ■ po (6(f)) 



(5a) 



where p + and po are the probabilities that, for a given node, the sum of its inputs 
is larger than or equal to the threshold 8, respectively. These probabilities can 



be written as (see Appendix B I 



P+ (b(t))= J2^n(h)J2 
k i = l i=0 



where L 



(1 - b) l b k 



E 

l=h 



(kj - i)a G + ' 
(a G + a R ) 



ki i 



p l q k ^- 1 , 



(5b) 



5 Given that statistical equivalence is assumed, then a will be representative of the 
entire network. Therefore, only the state of the inputs of <Ji is important, regardless 
of the states of all the other nodes. 



X 



/ ki\ ,~ V (^i A / ,/, . 



k i = l i=0 ^ 7 Z=0 ^ 

X ^a G i,a R (fe ! -i-0+f ■ ( 5C ) 

Since we are interested in the asymptotic value of the order parameter, it is 
necessary to compute the final number of active elements = lim^oo b(t). This 



is just the stable fixed point of the map given in Eq. (5a I, namely, = B (&oo)- 
From the set of Eqs. ([5]), the value of b^ is computed numerically for each 
particular realization of paramete rs. Once the value of boo is known, it is used 



to compute the influence l( kd K In Appendix C we show that l( kd ' is then given 
by 

K*d) - V" ( hl \ V" V f*\ ( h ~ l \ h l +™ n _ h \ k -- 1 - 



i=0 v 7 1=0 m=0 v 7 v 7 

x X(ki,kd,iJ,m), (6) 



where T{ki, kd, i, I, m) is defined as 



u f vf wf ( l\ (m\ ( A 

VtJ V v J V w ) \k d 



ki—i—m \ 
u—v—wJ 



( ki ) 

U=Uq V=Vq w=wq \kd' 

x {H (a G (l -u + w)- a R (m - v + z) - 9) ■ [(1 - b^ )5 aa i- aRmi g 
+H(anm + 9 - a G l)] + 5 aG ^_ u+w ^ aR ( m _ v+z)+9 ■ \h{t)5 aa i- aR m,9 
+ b OQ H{a R m + 9 — a G l) +(1 — b 00 )H{a G l — a R m — 9)] 
+ H (a R (m -v + z) +9 - a G (l - u + w)) ■ [boo8 aa i- aR m,e 
+H(a G l-a R m-9)}}. (7) 

where the summation is done between the limits of a multivariate hypergeometric 
distribution, [^] and H(x) is the Heaviside step function with H (0) = 0. 

Note that the influence l( kd > already contains information about the temporal 
correlations through the value of b^ . However, the above expressions are not 
exact because is computed from Eqs. [5j which were formulated using the 
mean-field assumptions. In spite of this approximation, it is an improvement 
over the original annealed approximation which completely neglects the temporal 
correlations. Once the value of I^ kd ^ if obtained from the above equations, it is 
used to obtain the Derrida map, which determines the temporal evolution of the 
Hamming distance, as 

00 ki / / \ 

h(t + 1) = M (hit)) = PM J2 I (kd) f! ) W)] kd [1 - Kt)p- k * . (8) 



6 Specifically we have that uq = max(0, kd + 1 — ki) , Uf = min(7, kd); vq = max(0, kd — 
u — (ki — I — m)), Vf — min(', kd — u); wo — max(0, kd — u — v — (ki — I — m — i + 1)) , 



wj = min(/, kd — u — v) 
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This equation tells us that the size of a perturbation avalanche after one time 
step depends on the probability to find k c i damaged input nodes between two 
configurations S t and S t , and on the probability I^ kd ^ that this damage spreads 
to the configurations S t +\ and St+i at the next time step. 



3.1 Sensitivity and influence of variables 



Once boo is obtained from Eq. (5a), it is possible to calculate the phase diagrams 
for the different parameters involved in a network realization using Eqs. ©, ^ 
and Q. There is however one last point that needs to be considered before the 
computation of the sensitivity, which is that the average influence of variables 
is not necessarily null. 

By definition, 1^ is the probability that, for a given node a„, a damage 
on none of its input elements spreads to the next time step. This means that 
<r n will have different values in the configurations St+i and E t +\ even when all 
of its inputs had the same values in the previous configurations E t and E t . In 
Kauffman RBN's this cannot happen, because the equality of the inputs of a n 
in the two configurations E t and E t guarantees that a n will have the same value 
in the next configurations E t +\ and St+i, and therefore, in this case 
However, for RTN's, the last line in Eq. ^ makes it possible for a n to be different 
in E t +\ and St+i even when all of its inputs were the same in the previous 
configurations E t and Et- This happens when \ a >n,j< J n :j (t) = 9 and a n had 
different values in the configurations E t and E t . In such a case a n will remain 
different in the next configurations St+i and E t +\. This can be considered as 
a damage spread for zero input variables, and consequently 1^ ^ 0. Note that 
this happens only when the equality in the last line in Eq. ^ is satisfied, which 
in turn occurs only for integer values ot 9. 

Taking the above considerations into account, it is possible to write as 

/ (0) = Po (boo) h(t) with 9 el, (9) 

Using the previous equation, one is finally able to get the sensitivity of the 
network, defined in Eq. Q, as 

oo 

S= PQ {boo)+ Kn{k t )hl {1) . (10) 

fei = l 

The derivation of the last two expressions is presented in |Appendix D| 

Note that 1^ also depends on fc^, and that both JW and po (boo) depend 
on the parameters of the network realization p, ag, olr and 9. Interestingly, the 
sensitivity S, and thus, the dynamical phase in which the system operates, only 
depends on the lower influences and I^. This means that the effect of small 
changes in the two configurations St and S t are the ones that determine the 



7 This is why in Ref. [29] the summation over kd excludes kd = 0, whereas in our 
Eq. Q the sum starts from kd = 0. 
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newtork dynamical regime. However, S t and E t cannot be arbitrary, since they 
must have a fraction of active elements close to the final one, b^. This restriction 
has profound effects on the initial apparent dynamical behavior of the network 
as compared to what actually happens at the end of the dynamics, in the sense 
that two trajectories that initially appear to converge may end up diverging, and 
vice versa. We will discuss this problem in Sec. [4] 



3.2 The phase diagram for the homogeneous random topology 



Eq. ( 10 ) determines the structure of the phase diagram as a function of the 
network topology (contained in Pi n (k)), and the other parameters of the system. 
Here we consider the homogeneous random topology Pi n {k) — 8x,k in which each 
node has exactly K regulators randomly chosen from anywhere in the system. 



In such a case, Eq. ( 10 ) establishes a relationship between the sensitivity S 
and the value of the parameters K, p, ag, and 9. The ordered phase occurs 
in those regions of the parameter space in which S < 1, whereas the chaotic 
phase occurs whenever S > 1. The critical region is the one for which S = 1. 
As the parameter space is 5-dimensional, an exhaustive exploration is neither 
illustrative nor computationally feasible. Instead, we present the phase diagram 
K vs. p for the following cases, which are representative of the general behavior 
observed across the entire parameter space: 

— Case 1: Activating and inhibiting interactions are of the same magnitude 
(a G = 1, a fl = 1); 

— Case 2: Inhibiting interactions are stronger than activating ones (og = 1, 
a R = 2); 

— Case 3: Activating interactions are stronger than inhibiting ones (og = 2, 
an. = i); 

— Case 4: Inhibiting interactions are always dominant (q,g = 1, clr — 20); 

— Case 5: Activating interactions are always dominant (<zg = 20, an = 1); 

Additionally, for each of the five cases listed above we used the threshold values 
6 = 0.5, 9 = and 9 = -0.5. 

The resulting phase diagrams are shown in Fig. [2] It is immediately apparent 
from this figure the asymmetric structure of the phase diagram with respect to 
the activator fraction (measured by p) and strength (measure by the quotient 
cig/cir). In general, it appears that activators strongly push the network into 
the frozen phase (in blue), while repressors move it towards the chaotic phase 
(in red), but less drastically. This can be seen in the extreme cases of dominant 
activators (ac = 20, an = 1) where the chaotic region almost dissappears, while 
for the opposite case of dominant repression (ac = 1, an. = —20) the frozen 
region is considerably smaller than the chaotic one. Another important point to 
note is the different behavior of the critical line for the three threshold values of 
interes: For 9 = 0.5 there are two critical values of p for each value of K, whereas 
for 9 — —0.5 and 9 = there is only one. In this sense, of all the cases shown in 
Fig. [2j the phase diagrams for 9 = 0.5 are the ones closer to the phase diagram 
obtained for RBN's [8]. 
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Fig. 2: Phase diagram p vs. K for threshold networks with different parameters, 



obtained by numerically solivng Eq. (10). The color code represents the value 



of the sensitivity S. The ordered phase (S < 1) is represented in blue while 
the chaotic phase (S > 1) is represented in red. Zones near the critical region 
(S = 1) are white, while the critical region itself is represented by the black line. 
Note the asymmetry of the phase diagrams, especially for 8 = and 9 = —0.5. 
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Finally, it is important to mention that we obtain the same results reported 
in Ref. [37] for the special case p = 0.5, ag = 1 and or = 1, but only for 
the threshold values 9 = 0.5 and 9 = —0.5. However, for 9 = we obtain a 
completely different behavior as the one reported in Ref. [27]. Indeed, we find 
that the phase transition occurs at K = 1, whereas the authors in Ref. [57] report 
that the phase transition ocurs between K = 12 and K = 13. This discrepancy 
is due to not properly taking into account the self-regulation conveyed in the last 
line of Eq. |2]), which happens only for integer threshold values (see Sec. 2.1 1. In 



the next section we present numerical results that support our analytic approach 
for integer threshold values. 

4 Numerical experiments 

To test the validity of the expressions obtained in Sec. [3] we performed numerical 
simulations of the network dynamics using ensembles of 100 RTN's, with TV = 
1000 nodes each. We used networks with homogeneous random topologies and 
varied K from K = 1 to K = 8. The other parameters p, ac Q>R and 9 were 
chosen to represent the different behaviors depicted in the phase diagrams shown 
in Fig. H 

4.1 The uncorrelated sensitivity So and the final avalanche size 

To compare the analytical expressions with the results of the numerical simula- 
tion we need to compute two parameters: The uncorrelated network sensitivity 
So and the final avalanche size . In Sec. 2.2 we describe how to compute 



for a given networkj:ealization. To compute Sq let us consider two initial con- 
figurations Sq and Eq that differ only in one element, namely, whose Hamming 
distance is 1/N: 

d(So,S ) = ^. (11) 



N 

Then, Sq/N is the average Hamming distance of the two configurations at the 
next time step: 

where the average (•) is taken over all possible pairs of initial conditions satisfying 



Eq. (11). In other words, So is the average size of the perturbation avalanche 
after one time step, given an initial perturbation of only one node. Note that 
Sq is the slope at the origin of the Derrida map without taking into account 
the correlations developed throughout time between network states. This is the 
reason why we call Sq the "uncorrelated" sensitivity. Therefore, for large N, Sq 



should be the sensitivity of the network given by Eq. ( 10 1 with &oo = £>o = 0-5, 
which assumes complete independence between the network states^] 



Since the initial configuration So in Eq. (Ill is chosen randomly from all possible 



configurations, the sequence of 0's and l's in Eq can be thought of as iV independent 
Bernoulli trials with probability 1/2, which gives 60 = 0.5 for the expected fraction 
of l's. 
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Fig. 3: Uncorrelated sensitivity So for RTN's with N = 1000 and different values 
of the parameters p, aa, a.R and 9. The symbols are numerical data computed 
using ensembles of 100 networks. For each of these networks, S was averaged 
over 10000 pairs of initial conditions differing in just one node. The error bars 
represent the standard deviation. The solid lines correspond to the theoretical 
result gien in Eq. (10) with b = 0.5. Note the excellent agreement between the 
theoretical prediction and the numerical data for all the different parameters 
used. 



Since the uncorrelated sensitivity Sq has no dependence on the correlations 
developed in time, we expect the analytic results derived in Sec. [3] to accurately 
reproduce the behavior of So. However, we do not expect this analytic compu- 
tation to describe as acurately the value of hoo , because in this computation 
the temporal correlations were approximately taken into account only through 
the value of , whereas the actual value of h^} depends on the precise way in 
which the network evolves in time. This is illustrated in Figs. [3] and [4j 



It can be seen in Fig. [3] that the uncorrelated sensitivity Sq computed numer- 
ically (symbols) shows a remarkable agreement with the theoretical prediction 
(solid line) for the different combinations of parameters used. It is interesting to 
note the variety of behaviours exhibitted by So in RTN's, which is in marked 
constrast with the linear behavior observed in standard Kauffman Nets. Indeed, 
for RBN's Sq — 2p(l — p)K, whereas for Kauffman networks with canalyzing 
functions S = l/2 + (if-l)/4 [29 . In contrast, Fig.[3]shows that for RTN's the 
dependance of So on K is nonlinear and can even change inflexion or decrease 
with increasing K. This general nonlinear behaviour occurs even for the simple 



cases p = 0.5, ac = Or = 1, = 0, and 9 = ±0.5, where we have (see Appendix 
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[E| for a derivation) 



K + 1 (2K 



2 2K \ K J im9 = 

So={ ,\ (12) 

in which So ~ V K for large if. The above result allows the network to ramain 
close to the critical phase for a wider range of values of K than standard RBN's. 
This might be important given that there is evidence showing that real genetic 
networks, in which the gene input connectivity varies considerably from one gene 
to another, work near the critical phase Sq = 1 [34] . 

With regard to the final size of the perturbation avalanche, Fig. [4] shows the 
value computed numerically (line with symbols) , and the value h*^ predicted 
by the analytic computation of Sec. [3] (solid line). Although and h*^ are 
qualitatively very similar to each other for the cases 9 = ±0.5 depicted in Fig. |4j 
their quantitative correspondence is not as precise as it was for Sq. As it was 
mentioned before, this lack of precision was expected due to the approximation 
in the computation of the temporal correlations. However, for integer threshold 
values /loo and h*^ do not necessarily agree even qualitatively, as it is illustrated 
in Fig. [l] for the special case 9 = 0. We discuss the origing of this discrepancy 
further below. In the mean time, it is important to emphasize the reason why 
the Derrida map is not always useful to discriminate the dynamical regime in 
which the network operates. 

Fig. [6] shows the temporal evolution of the average Hamming distance h(t) 
between two trajectories that started from two slightly different initial conditions 
Sq and Sq. (The average is taken over many pairs of initial conditions.) In all 
the cases shown in this figure, h(t) decreases in the first time steps. Therefore, 
according to the Derrida map, which takes into account only the first time step, 
these networks should be in the ordered regime. However, after that initial de- 
crease, the Hamming distance h(t) increases again reaching a value considerably 
larger than the initial Hamming distance h(0). Thus, in the long term the initial 
perturbation is amplified, which means that the dynamical regime in which these 
networks operate turns out to be chaotic. It should be noted that the behavior 
reported in Fig. [6] was obtained for networks with biologically reasonable values 
of the parameters: clq = clr = 1, 9 = 0.5, p = 0.5 and K = 3. 



4.2 The = case 

We now address the anomalous case 9 = 0. As we have seen in the previous 
section, in this case the annealed approximation gives very accurate results for 
the initial sensitivity Sq but very poor results for the final avalanche size h& ■ 
As discussed in Sec. [2j integer thresholds allow the possibility for some nodes 
to become frozen whenever their input sum in Eq. ^ equals the threshold. 
These frozen nodes generate explicit correlations in time, which in turn produce 
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Fig. 4: Final size of the perturbation avalanche computed numerically (hoo , red 
circles), and analytically as the fixed point of Eq. ^ (K>oi sou d black line). The 
numerical data were computed for an ensemble of 100 RTN's with N = 1000 
and different values of the parameters p, a,Q, clr and 9. The important point in 
this figure is the use of noninteger threshold values: 6 — ±0.5. For each network 
realization, we averaged over 10000 pairs of random initial conditions with 
Hamming distance ho = 0.1. Note the large standard deviations in the numerical 
dada (error bars) . Despite this enormous variability in each network realization, 
the average numerical data qualitatively follow well the theoretical prediction. 



a strong dependance on the history of the dynamics, and thus, on the initial 
conditions. This is illustrated in Fig. [fj] where the final avalance size h^} is 
plotted against the initial perturbation size h = h(0) for networks with p = 0.5, 
a G — a R = 1 an d = 0. It is apparent from this figure that for K < 5 the value 
of h£a strongly depends on ho, and this dependece becomes less strong as K 
increases. This is because for large values of K it is harder for the input sum to 
equal the threshold. 

Another problematic consequence of using integer threshold values is the 
existence of an enormous number of punctual attractors, many of which differ 
only in the value of just one node. This anomaly has been noted before in 
Ref . [27] . It also occurs in the cell cycle models of S. cerevisiae [15] and S. pombe 
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Fig. 5: Temporal evolution of the Hamming distance for 5 different random 
threshold network realizations with N — 1000, ao = &fl = 1, = 0.5, p = 0.5 
and K = 3. Each curve is the average over 10000 randomly chosen pairs of ini- 
tial conditions separated by a Hamming distance ho = 0.1. Note that initially 
the Hamming distance decreases. Using only the Derrida map, which takes into 
account only the first time step, one would conclude that the networks operate 
in the ordered regime. However, the correlations developed in time due to the 
network structure and the number of active nodes make the Hamming distance 
rise again and approach a nonzero value, which is characteristic of the chaotic 
regime. 




Fig. 6: Final size of the perturbation avalanche as a function of the initial 
perturbation size tiQ. Each point is the average over 100 random threshold net- 
work realizations with N = 1000, aa = clr = 1, p = 0.5, 6 = 0, and 1000 pairs of 
random initial conditions for every ft, i n each of these networks. Note the strong 
dependence of n£) on ho, especially for small values of K where the temporal 
correlations are stronger. For large values of K these correlations become weaker 
and consequently h&) becomes almost independent of h . 
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Fig. 7: Average number of attractors as a function of the network connectivity 
K for RTN's with (a) N = 100 and (b) N = 15 nodes. In each panel, we 
used p — 0.5, ac = a.R = 1 and 6> = (circles), = 0.5 (squares) and 9 = —0.5 
(triangles). For the large networks in (a) we sample the configuration space using 
50, 000 randomly chosen initial conditions, whereas in (b) the full configuration 
space was probed. Both in (a) and (b) each point is the average over ensembles of 
50 networks. Note the extremely large number of attractors obtained for 8 = 0, 
especially for moderately small values of K. In particular, for 9 = in (a) almost 
every sampled initial condition leads to a different attractor. 
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Fig. 8: Structure of the attractor landscape of RTN's with N = 15, K = 8, 
p = 0.5, a_R = flg = 1, and 9 = 0.5 (left) and 9 = (right). Each of the basins 
of attraction shown is the largest one in a given network realization. Note that 
for 9 = 0.5 the attractors have several configurations (black dots at the center 
of each structure), whereas for 9 = all the atractors are punctual (only one 
dot at the center). Additionally, for 9 = 0.5 the long arm-like structures indicate 
the existence of long transients and a reduced number of configurations which 
all the routes that go to the attractor have to go through. Contrary to this, for 
9 = the basins of attraction look more sparse and with shorter transients and 
more distributed across the configuration space. 
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|16j , where many punctual attractors with small basins of attraction were foundj^j 
Fig. [7] shows the average number of attractors as a function of the network 
connectivity K for 8 = (circles), 8 — 0.5 (squares) and 8 = —0.5 (triangles). 
In all cases we used p = 0.5 and ac — a,R — 1< Fig. [7£i corresponds to large 
networks with N = 1000. In this case the number of possible configurations 
is astronomically huge (Q = 2 1000 ). Therefore, an under sampling of the state 
space has to be done, in which case only some attractores will be found. We 
sampled 5 x 10 4 configurations. Surprisingly for 8 = almost every sampled 
initial configuration ended up in a different attractor. The same happens for 
smaller networks with N = 100, as it is shown in Fig. [TJd. However, the same 
undersampling performed in networks with non-integer threshold values (8 = 
±0.5) reveals a number of atractors which is several orders of magnitude smaller 
than the one obtained for 8 = 0. Finally, Fig. |7J; shows similar results but for 
small networks with N = 15, for which the entire state space can be probed 
(J? = 2 15 = 32768). Note that in this case, the average number of attractors 
decreases with K for 8 = 0. This behavior is marked contrast with the one 
observed for non- integer threshold values (and for RBN's), where the average 
number of attractors grows with the network connectivity K. 

The 8 = case presents "anomalous" behavior not only with regard to 
the number of attractors, but also in the structure of the state spacej^] Fig. [8] 
shows the largest basin of attraction for three network realizations with ./V = 15, 
K = 8, p = 0.5, an = aa = 1, and 8 = 0.5 (left part) and 8 = (right part). For 
those parameters the networks are in the chaotic regime. It can be seen from 
this figure that the structure of the largest basin of attraction for the non-integer 
threshold is characterized by long transients and attractor length. This is similar 
to what is obtained using RBN's with the same p and the same K . However, 
for 8 = the structure is quite different. Note first that all the attractors are 
punctual. Although this is not the rule, it is the most probable situation for 
8 = 0. Additionally, the transients are comparatively short and the whole basins 
look somehow sparse as compared to the ones on the left part. The long arm- 
like structures observed in the basins of attraction for 8 = 0.5 reflect that the 
routes to reach the attractor are concentrated in a few number of states. From a 
biological point of view, these few states can be considered as the "checkpoints" 
of the differentiation (or metabolic) pathway. Contrary to this, the sparseness 
observed in the basins of attraction for 8 = indicate that the routes to reach 
the attractor are much more distributed across the state space, and therefore, 
the existence of "checkpoints" is harder to attain. 



However, in these cell cycle models the authors deemed the attractors with small 
basins of attractions as biologically irrelevant and neglected them. 
Here, by "anomalous" we mean with respect to what is observed in standard Kauff- 
man networks. 
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5 Summary and discussion 

We have investigated the dynamical properties of random threshold networks 
(RTN's), which differ from standard Kauffman networks (or random Boolean 
networks, RBN's) in that the regulation of the state of the nodes is done by 
means of threshold functions. These networks have been used in the modeling of 
genetic regulatory networks of real organisms using parameter values that seem 
biologially meaningful |15I16| . An important characteristic of the threshold net- 
work model is that it assumes that gene regulation is an additive process. Namely, 
that the combined effect of the regulators of a given gene on the state of that 
gene is just the sum of the positive regulations minus the negative ones. Because 
of its simplicity, this assumption is very tempting when constructing models 
of gene regulatory networks. However, there are many examples in real sys- 
tems showing that gene regulation is combinatorial rather than additive, which 
means that the effect of some regulators (i.e. whether activatory or inhibitory) 
depends on the presence or absence of some other regulators^] In several cases, 
these combinatorial processes are represented by Boolean functions that cannot 
be obtained from threshold functions (for instance, the XOR function), which 
makes the additivity assumption mentioned above questionable. 

Additionally, a deeper analysis of the dynamics of RTN's reveals anomalies 
inconsistent with the expected behavior of gene regulation models, precisely for 
the "biologically meaningful" values of the parameters that have been used. Of 
particular importance is the case of integer threshold values (illustrated here us- 
ing 8 = 0), where the networks typically have an enormous amount of attractors. 
Also in this case, there is a sharp disagreement between the ensemble properties 
predicted by the annealed approximation and the ones observed in the numeri- 
cal simulation using concrete network realizations. It is worth emphasizing that 
this disagreement happens neither for RTN's with non-integer threshold values 
nor for RBN's. It is not surprising to find such a disagreement in specific net- 
works constructed in very peculiar ways (as the ones used for the modeling of 
the cell-cycles). What is surprising is that the typical members of the ensemble, 
constructed in a completely random way, present such anomalies. In fact, the 
cell-cycle networks in Refs. |15ll6j do exhibit these anomalies, as they have a 
large number of attractors. However, the authors of that work considered most of 
these attractors as biologically irrelevant because of their small basins sizes, and 
neglected them. Nonetheless, from an evolutionary point of view it is not clear 
whether or not the size of the basin of attraction is relevant, as it is not known 
whether this parameter is under selective pressure. Actually, in other studies of 
gene regulatory networks of real organisms, the biologically meaningful attrac- 
tors of the wild-type organism do not possess the largest basins of attraction, 
but on the contrary, sometimes they have very tiny basins [12113] . 

We computed analytically the phase diagram that determines in which pa- 
rameter region the network has chaotic, ordered or critical dynamics. Contraty 
to what happens for RBN's, the phase diagrams obtained for RTN's are always 



11 An common example of this are dual regulators in E. coli [41 j . 
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asymmetric with respect to the fraction of positive regulations p. It is only for 
non-integer positive thresholds (illustrated here for the case 8 = 0.5), that the 
phase diagram looks semi-symmetric, similar to the one obtained for RBN's (see 
Fig. [2| . This may be important because in such a case there is a bigger freedom 
to vary p and still remain close to the critical region, especially at low newtork 
connectivities as the ones reported for networks of real organisms that operate 
in the critical regime (A" ~ 2, [34 ). Note that the case 9 = 0.5 is biologically 
reasonable not only in terms of the phase diagram, but also because it corre- 
sponds to a situation in which at least one positive regulator has to be active in 
order to activate the target gene. Contraty to this, in the case 9 — —0.5 all genes 
get activated when all of its positive become innactive, which is an artifact of 
the model rather than a biological observed behavior. Furthermore, the phase 
transition for 9 = —0.5 and p — 0.5 occurs at a network connectivity K = 4 (seel 
Fig. [4]), which is large compared to the one observed in real networks. Therefore, 
this case with negative threshold values seems to be inadequate for the study 
of the theoretical properties of gene regulatory networks. Consequently, some 
of the conclusions about the evolution of RTN's with negative thresholds might 
have to be reinterpreted [43] . 

One important point analyzed in this work was the usefulness of the Der- 
rida map to elucidate the network's dynamical regime. As it is shown in Fig. [6j 
for RTN's the first steps in the dynamics may indicate that the network is in 
the ordered regime, while in fact the long-term behavior is chaotic. This oc- 
curs when long-term correlations are developed thoughout time, which always 
happens in RTN's, especially for integer threshold values. For in such a case, 
the self-regulation implied by the last line in Eq. ^ induces long-term mem- 
ory in the system. For RBN's these long-term correlations do not exist, and 
the Derrida map accurately predicts the network's dynamical regime. This is 
an important aspect that has not been properly taken into account in current 
work that aims to characterize the network's dynamical regime in real biological 
networks. Therefore, a more thorough study is necessary in this direction. 

Finally, it is important to note that in this work we used the same fixed 
connectivities, interaction strengths and thresholds for all genes. However, more 
realistic situations would require assigning these quantities differently to the 
different genes in the network. For instance, for some genes the inhibitory regu- 
lators may be dominant, whereas for other genes the activatory regulators would 
dominate. Also, genes with integer as well as non-integer threshold values can 
coexist in the same network. Exploration of these possibilities can reveal dynam- 
ical behaviors more consistent with biological systems, which in turn would help 
to discern the model's characteristics relevant for biological modeling. 
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Appendix 

Appendix A: Derivation of the map b(t + 1) = B (b(t)) 

We first remember that b(t) represents the fraction of active nodes (<7j = 1) for 
a given network configuration at time t. Using the annealed approximation, the 
evolution of b(t) depends only on the previous state and can thus be given by 
the map B, which relates the number of active nodes after two consecutive time 
steps. 

To explicitly obtain B let us consider the following. Since in the annealed 
approximation we assume statistical independence between the nodes, the frac- 
tion of active elements b(t) can actually be considered as the probability that 
an arbitrary node cr, is active at time t. Therefore, b(t + 1) corresponds to the 
probability that a node is active after one time step. From the dynamical equa- 
tion for the nodes, Eq. ([2]), it is apparent that there are only 2 ways in which 
this can happen: Either the sum y\ a n ja nj (t) was larger than 9 or it was equal 
to 9. In this last case we additionally need the node itself to be active at time t 
so that it is still active at time t + 1, which happens with probability b(t). If we 
denote p + {b(t)) as the probability that yj ■ a n,iV nj {t) > 9 and po (b(t)) as the 
probability that y_^ . a n jO nj (t) = 9, then B must be the sum of the probabilities 
of these two events: 

b(t + 1) = B (6(t)) = P+ (6(t)) + b(t) ■ po (b(t)) , (1) 

This corresponds to Eq. ( |5a[ ) in the main text. The full expression for po and p+ 
are derived in the next section, |Appendix B| 



Appendix B: Derivation of po (&) and p+ (b) 

To derive these expressions we use the mean-field method from Ref. [53]. Let us 
denote Ps(y) as the probability distribution function of the sum = Y^j=i a ij a H 
in Eq. ^ of a node Oi with ki inputs. The probability that £j = and £j > 
are then given, respectively, by 

p + (b,ki) = ]im P s (y)dy (1) 
p (b, ki) = lim/ Ps(y)dy. (2) 

The weights a^j can be considered as random variables which take the value ac 
with probability p and —an with probability q = 1 —p, as defined in Sec(2] Using 
this and denoting b as the probability that a node is active, we can consider the 
£y = dijOi. as random variables which can take the values with probability 
1 — b, a<3 with probability bp and —an with probability bq, that is 



P s (x) = (1 - b)S(x) + bp S(x - a G ) + bq 5(x + a R ), 



(3) 
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Using the statistical independence assumption of the annealed approximation, 
each £jj in the sum £j = X^jLi £ij i s an independent random variable with 
probability distribution P^(x). Because of this, is the sum of fc f independent 
random variables, and thus Ps{y) must be the fcj-fold convolution of P^: 



P E {y)=P t *P t *...*P t {y). 



ki times 



Taking the Fourier transform of the above equation we get 



(4) 



Pz{co)= P^cu) 



(5) 



where P e = (1 - b) + bp e - lu)aa + bqe iujaR . Thus, P s (y) is obtained by taking the 
inverse transform of the last equation. Using the binomial theorem twice we get 



ki 



2TT J_, 



2ir 
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(6) 



X / e «^y e -iua G l e iua R (ki-i-l) ^ 
J — oo 

ki ki — i 

( 

i 



= E E ( ^ ) H byb k ^ P l q k ^- l 8 [acl - ante — i — I) — y] 

(7) 
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If we substitute this last result into Eqs. ([!]) and ([2| we find 

/•OO 

P+ Q>,k) = lim / Ps(y) dy 



i=0 x 

X 



lim£° | £ J ^ (?*•-*- <<5 [a G Z - a^fe - i - l) rfy 

= e ft) a - e ft r 4 ® 

i=0 ^ ' l=U ^ 7 

where u = {k *~ l)aG + e + i 

a G + a R 



p (b, h) = lim / Ps{y) dy 



'k, 



i=0 



E(7J(i-^- 

9+e f k i~ 



x lim jf j E ; ^g**-*-'* M ~ - * - - y} \ dy 

- e ft) a - e ft r Vfl fc, - < - , *«i,.-(*.-i-o + « ( 9 ) 



i=0 v 7 Z=0 



where in Eq. pj the minimum value of Z = ij was chosen so that the argument of 
the Dirac delta function is always above 9, as specified by the limit. Similarly in 
Eq. ^ it is chosen so that it is exactly equal to 9. Finally, since the probability 
distribution of fej is given by Pi n (k) we have that 



p+ (b, hi) = ^ Pin(h) P+ (b, h) 



fci=i 

oo hi /, v ki — 



E E ( 1 ) a - E ( fcl ; >V i_i ~'' ( 10 ) 



fc 4 =l i=0 v 7 1=1 

oo 



Po 



= ^2 Pin(ki) po {b,ki 



k i = l 

oo ki /, \ ki — i /j 



plqkt-i-l 



fc 4 =l i=0 ^ 7 2=0 
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which correspond, respectively, to Eqs. ( 5b I and (5c I. 



Appendix C: Derivation of for Boolean threshold networks 

We first remember the definition of I^ kd \ the influence of kd variables. Denoting 
S t and S t as two configurations in which kd of the inputs of an arbitrary node 
<Xj are different, I^ kd ^ is the probability that, after a time step, the node <Ji in 
the new configurations, S t +i and S t+ i, are different from each other. 

Consider the average over all possible active/inactive and activatory/inhibitory 
configurations of the inputs of an arbitrary node <7j. If the node has fcj inputs, 
activator probability p and a fraction of active nodes b, then this average is given 
by 

- E ('/)"'«" 'EE (T) (! - bf>- l - m X. (1) 

i=0 ^ ' ;=o m=0 V / \ / 

Here p l q ki ~ l (g = 1 — p) is the probability that for a given input configuration 
with fci regulators there are i activatory interactions and ki — i inhibitory ones, 
which can be chosen in ( .*) possible ways. 

b l+rn ^ _ tyki-l-m j g probability 

that there are I active activatory inputs and m active inhibitory ones, which 
can be arranged in ( ki r ^ t ) different ways. Since I^ kd ^ is the probability for one 
arbitrary node (regardless of the number of active or inactive inputs), we have 
to compute the average over all possible input configurations of Tiki, kd, i, I, m), 
which is the probability that a damage spreads when kd of the input elements are 
damaged given that this configuration has i activatory inputs and fcj — i inhibitory 
ones, which in turn have I and m active/inactive input nodes, respectively. 

To find X(ki,kd,i,l,m) we need to consider all possible ways in which the 
damaged nodes can be arranged. There are / active activatory input nodes and m 
active inhibitory ones, and thus, i — I inactive activatory elements and ki — i — m 
inactive inhibitory ones. Therefore, we may have u damaged active activatory 
inputs, v damaged active inhibitory ones, w damage inactive activatory ones and 
z = k d -u-v-w damaged inactive ones. Since there are ( l )( m ) (, k *- l ~™ ) 

possible ways in which damage can be distributed, and given that there are (^*) 
total ways in which the damaged nodes can be arranged, then the probability 
for each value of it, v, w is given by a multivariate hypergeometric distribution 



I \ ( m\ I i — l\ ( ki — i — 



, sitjyvjy w J\kd-u-v-wj 
Pt(u,v,w) = y—^ v = Vo,...,Vf , (2) 

Ki 



k d 



U = U ,...,Uf 
V = Vo,.-.,Vf 
W = W ,...,Wf 



where 



uq = max(0, k d + I — ki) , Uf = min(7, k d ) 

Vq = max(0, kd~ U — {ki — I — m)), Vf = min(7, kd — u) 

w = max(0, kd — u — v — (ki — I — m — i + I)), wt = min(Z, k d — u — v). 



(3) 
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Finally, we need to consider all possible ways in which a damage can actually 
make the state of <7j in the damage and the undamaged configuration different at 
time t+1. From the definitions of I and m it follows that, before damage, there 
are I active activatory input nodes and m active inhibitory ones. After damage, 
using the definitions of u, v, w and z, there will be I — u + w active activatory 
input elements and m — v + z active inhibitory ones. Using this information and 
Eq. pi), it is clear that the damage can spread in 3 different ways: 



i) If in the damaged configuration the sum is above the threshold (ac(/ — u + 
w) — a R (m — v + z) > 9), then the state of the nodes can be different if, 
without the damage the sum is either below the threshold {a G l — a R m < 9) 
or exactly at the threshold {clqI — a R m — 9) but with the condition that it 
was inactive at the time-step before (so that they are different in the next 
time step), which happens with probability 1 — 6. 

P\ = H (a G {l -u + w) - a R (m - v + z) - 6) ■ [(1 - b)5 aG i- aRm fi 

+H(a R m + 9- a G l)} . (4a) 



ii) The second possiblity is that in the damaged configuration the sum is just at 
the threshold (ag(( — u + w) — a R {m — v + z) = 9). In that case damage can 
spread if: (a) before the damage the sum is above the threshold (acl — a-Rjn > 
9) but only if it was inactive before (with probability 1 — 6); (b) if before the 
damage the sum is below the threshold (acl — clrtu < 9) but only if it was 
active on the step before (with probability 6) ; and (c) if before the damage 
the sum is again exactly at the threshold (clqI — a R m = 9) but only if both 
configurations where damaged before (with probability h(t)). 

P-2 = 5a G {i-u+w),a R {m-v+z)+e [h{t)5 aG i- aRm fi + 6 H(a R m + 9 — acl) 
+ (1 - b)H(a G l - a R m - 9)] . (4b) 



iii) In the third possibility the damaged configuration has its sum below the 
threshold (aa(l— u+w)— a R (m— v+z) < 9) so the state of the nodes will differ 
if before the damage the sum is either above the threshold (agl — a R m > 9) 
or exactly at the threshold (clqI — a R m = 9) but only if it was active at the 
time-step before, which happens with probability 6. 

P 3 = H {a R (m -v + z) + 9 - a G (l -u + w))-[b 5 aG i_ aRmfi 

+H(a G l - a R m - 9)} . (4c) 



Since all three cases can make damage spread, then the total probability for 
damage spreading X is the sum of these three possibilities: Eqs. (4a), (4b) and 
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(4c), averaged over all possible arrangements of damaged configurations, Eq. ([3]), 

Uf Vf Wf 

1= E E E Pr(u,v,w)-(P 3 + P 2 + P 3 ) 



U—Uq V — Vq W — Wq 

u f v f w f /l\/m\/i-l\/ ki-i-m \ 
\u) V v / V w J Kk^ — u—v — wJ 



EE E 



U—Uq V—Vq W—Wq V^d' 

x {H (a G (l -u + w)- a R (m - v + z) - 9) ■ [(1 - b^ )5 aG i- aRmi g 
+H(a R m + 9 - a G l)} + 5 aG {i-u+w),a R (m-v+z)+6 • [h(t)5 aG i-a R m,8 
+ b OQ H(a R m + 9 - a G l) +(1 - b 00 )H(a G l - a R m - 9)] 
+ H (a R (m -v + z)+9 - a G (l -u + w)) - [&<xA o J-<wn,0 
+H(a G l-a R m-9)}}. (5) 

Finally, averaging this damage spread for a given input configuration kd,i,l, m) 
over all possible input configurations using Eq. 0, we get the average influence 
of fcrf variables: 



I (k ^ = (l(ki,k d ,i,l,m) ) 



ic 



t (";)'' '< 'EE fi)f A '„ i y +m o >>■ ^ 



, i J — ' — ' \l I \ m 

i=0 v ' 1=0 m=0 



xX(ki,k d ,i,l,m). (6) 

Eqs. ^ and ^ with b — b^ correspond to the formulas and ([6]), respectively, 
in the main text. 

Appendix D: Derivation of S and JW 



As discussed in Sec. 4.2 is the probability that, for an arbitrary node <7j, 
damage spreads at the next time step when none of its input elements are dif- 
ferent between the two initial configurations, U t ,S t . Because of the possibility 
of having this sum giving exactly the threshold value 6>, is zero only for 
noninteger thresholds. 

The case for integer values of 9 can be obtained from Eqs. ^ and 
However, from Eq. ([2| we can see that the only way for a damage to spread 
when none of the input elements are different is by having J2j a n.jO'n j (t) = @- 
Thus, must correspond to the probability po that the sum gives exactly 



the threshold value, Eq. (5c). Since we need to consider the nonergodicity of 
the system, we use the final fraction of activatory nodes as the value of b. 
Now, for damage to spread not only does the sum need to be at the threshold, 
but also the two nodes must be different initially, otherwise they would be the 
same in the next time-step and the damage would not spread. Since h(t) can be 
considered as the probability that two arbitrary nodes are different, then 
must be the multiplication of both probabilities 

/ (0) = M&°°) • h(t) withfleZ, (l) 
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Taking this last result into consideration we can now calculate the sensitivity S 
of RTN's. Using g in @ we get 



S — ^ ' Pin{ki) 



ki=l 



dh 



+ fc .(jCi)_ 7 (0)) 



(2) 



h=0 



From Eq. Q it is apparent that I^ 1 does not actually depend on h. This happens 
because the middle term of Eq. ^ , which is the only one with a h dependance, 
can never be nonzero for the case of — 1. This is because refers to 
the case where after kd changes of an arbitrary input, the summation of both 
configurations is still at the threshold, which cannot happen when only 1 input 
is changed. With this information and using Eq. in ([2| we finally find 



S = P ^( k i) [Po ( b oo) + h - I {0) ) po (boo) • h 

kt = l 
oo 

= Pin(ki) (p (b QO ) + k t I ( - 1 
fci=i 

oo 

= P0 (boo)+Y, Pin{h)kil (1) . 



h=0 



(3) 



ki = l 



Eqs. ([T]) and ^ correspond, respectively, to Eqs. ^ and ( 10 ) in the main text. 



Appendix E: Derivation of for p — 0.5, aa — cir — 1, 6 — and 
e - ±0.5 



We first remember from Sec. |4.1| that So, the uncorrelated network sensitivity, 
is the average number of nodes by which two configurations differ after one time 
step if they initially differed in only one element: 



S = N(d(E 1 ,S 1 )), with d(E Q ,£ ) = 



1 

N' 



Because of the thermodynamic limit assumed in the annealed approximation, 



So should correspond to the sensitivity of the network given in Eq. ( 10 1 with 
&oo = bo = 0.5. This last choice of b is a consequence of the initial configurations 
being chosen randomly. In what follows we consider consider the case ni which 
Pin(k) = Sk,K with p = 0.5, a G = a R = 1 } 6 = and 9 = ±0.5. 



E.l S for p — 0.5, a G — a R — 1 and = Using Eqs. Q and |To| for 

the integer threshold 9 = 0, we have 



So=Po (b ) + KI^. 



(1) 
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where from Eqs. |7| and (5c 



1 



j=0 V / Z=0 

i K-i 



K-i 
I 



>l,K-i-l, 



i=0 v 7 2=0 m=0 v 7 v 



(2) 

= 0,6 = 0.5). 

(3) 



For T{o,q = 1, or = 1, = 0, b = 0.5) we consider the following. Since we have 
an integer threshold 8 = and we are looking for the influence of 1 variable, the 
input elements of the damaged and the undamaged configurations differ only 
by one node. Additionally, since the positive and negative weights have equal 
strenght, the only possible way for the sum to change sign is if: (a) without 
damage the node is exactly at the threshold (the Kronecker delta terms in Eqs. 



4a and Eq. 4c I , or (b) if after the damage it its exactly at the threshold (the first 



Kronecker Delta term in Eq. 4b ) . Since by damaging a single node we are not 
able to have the node again at the threshold because of the weights, the term 
with h(t) in Eq. 4b cannot be attained. Finally, since kd = 1, the only possible 
values for u, v, w and z = k ( i — u — v — w are 1 on one of them and on the rest. 
Using all this information and Eq. [5] wc have 

1 



1 = 



Uf Vf 



Wf 



2 ^ 

U = Uq V=Vq W=Wq 



Pr(u, v, w) {H{1 — u + w — m + v — z)5i t , 



+5i-u+w, m -v+z [H(m - I) + H(l - m)] + H(m - v + z- l + u~ w)8i im } 



m 
K 



I 



+5l-l,m 
1 



I 

K 



K 

K — i — m 



Si 





m 


i-V 




-1 


K 


+ K _ 




" I 


_ _L_ 


K-i- 


m 


K 


K 





2K 



K 

[S m ,i-i (K-i + l) + 5 m . l+1 (i + 1)] + ^Si 



(4) 



E.l.l po (bo) To get the result from Eq. ([2| we use the finite Laplace Trans- 
form method. Let us define a generating function 



fe=0 



k\ 



with 



i=0 



c k = 2 2k Po(MlK= fc = EL-P ^ fc 

k—i 

/(fc-o=E 



1=0 



k — i 
I 



Si.k- 



k—i—l' 



(5) 

(6) 
(7) 
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Using Eq. ^ and the definition of given in Eq. ^ , we obtain 



Aj=0 
oo k 



fc=0 i=0 



With the change of variable u — k — i and inverting the order of summation we 
get 



oo oo 



EE^-'VW 

j=0 tt=0 



E 

i=0 



i! 



E^/c 



Substituting into this last result the value of f(u) given in Eq. ([7]), making the 
change of variable v = u — I, and again exchanging the sums we get 



g(z) = e 2 



.tt=0 1=0 y ' 



oo_ oc_ 



EE 

i=o v=a 
e 2z Io(2z), 



IW. 



(8) 



where Iq is the modified Bessel function of the first kind and where we used the 
series representation of the exponential function and Io(z) = TJ°^ (z/2) 2t /(i!) 2 . 
Using the integral representation of Iq(z): 



Io(z) 



2tt 



and the change of variable 4> = 28 in Eq. (18]) we get 



1 

2n 
1 

47T 

oo 

E 

fc=0 



3 2z(cos8+l) dg 



— TV 

2tt 



O 4zcos 4> 



2tt 



IT 



k f27T 



4tt 



„2fc 



2tt 



(9) 
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Comparing this last result with Eq. (|5| and using the formula 



71/2 2k j. j j. (2fc-l)l! 1 {2k 



(2fc)H 2 2fc Vfc/' 



we finally find 



4 fc 


,.2-k 


iw 


-2tt 








JO 


/2k^ 




\ k , 


)■ 



c k = -r- I cos 2k ( 



cos 2fc ( 



which, using Eq. gives us the first term of Eq. ([!]): 



1 /2X 



(10) 



Po(bo) = W [ K )- (11) 



E.1.2 From Eqs. Q and @ we define 



= ~ (A + / 2 + /a) , (12) 



where 



i=0 v 7 1=0 m=0 v 7 v 7 

a- it (f) EE QC" ! )w*+i). <», 

i=0 v 7 (=0 m=0 v 7 v 7 
i=0 v 7 1=0 m=0 v 7 v 7 

To reduce these expressions we will use Vandermonde's identity and the mean 
value of the hypergeometric function 

m+ " A -t (:)(",) <«> 



r 



. fc / \r — k 

k=0 v 7 v 



fc=0 



Vfc("M( »M = jW'» + 'M (17) 

^ — ' \ b \ r — hi m 4- n \ r / v/ 



k J \r — k J m + n 
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Using these expressions in /i with the change of variables I' = 1 — 1 and i' = i — 1 

K 



i=Q 
K 



1 

2K 



i=0 
K 



(K-i + l) £ 

1=0 
i-1 



i\fK-i 
i 1-1 



E 7 c*-' +1 > £ 



1 , , /if 



2if 



i=0 



i'=0 

i - 1 



if -i 
V 



»'=0 



if — 1 — i' 







^2if\ 




-i 


v ^ / 



(18) 



Doing something similar for ji with the change of variable ml = m — 1 



2A" 



i=0 
K 



eO+d e 



if-? 



2if 



i=0 
K 



m=0 
Jf-t-1 



m — 1 



K - i 
rn 



m'=0 



K-i 

m'J\K — i— 1 — m' 



i=Q 



2if 

i=0 

if + 1 / 2if 



if 

i 

if 



if 

K-l-i 
if 

K-l-i 



4if Vif-1 



1 /2if 
4 V if 



(19) 



For /3 we just need to use Eq. ( 16 1 with m = n = K 

/ i\ I if — i 



1 



A" 



i=0 



Elf) E 



= ^E(f)E(,. 

2=0 V 



i=0 

1 K 

-Y 

2 ^ \ i 

i=0 



i=0 
A \ /if 

i 



K\ ( if 



if -i 



1 /2if 



if- j/ 2 V if 



(20) 



XXXV 



Using Eqs. Q, (19} and ([20} in §V2$ we find the second part of Eq. 0: 



I« = 



1 /2K 
2™\K 



(21) 



We finally get the first part of Eq. ( 12 1 using Eqs. (Ill and ( 21 ) in Eq. ([!]), which 
gives 

K + 1 (2K> 



S (p = 0.5, a G = 1, a R = 1,0 = 0) = 



2 2 ^ V # 



(22) 



E.2 So for p — 0.5, ac = = 1 and 6 — 0.5 In this case the threshold 
is a noninteger value, 9 = 0.5. Therefore, from Eq. ^ and (10), we have 



So = KI^. 



(23) 



where from Eq. (5c 



I' 1 



1 

%2K 



K 

E 

i=0 



K 



i K-i 

EE 

2=0 m=0 



K 



l ^l(a G =l,a R = l,6 = 0.5, b = 0.5). 

(24) 

In order to calculate I(ac = l,dR = 1,6 = 0.5, b = 0.5) we consider the follow- 
ing. The noninteger threshold 9 = 0.5 makes all the terms with Kronecker deltas 
effectively zeros, since the exact value of the threshold cannot be attained. Given 
that 9 = 0.5 > 0, it also makes nodes in which the sum of the updating func- 
tion equals 0, ■ a n ja nj = 0, become inactive, as this sum is smaller than the 
threshold. As a consequence, the only way in which the changing of one of the 
inputs of a node changes the state of the target node is if either the sum before 
the damage was and after the damage it is above 0, or vice versa. In addition, 
since kd = 1, then the only possible values for u, v, w and z = kd — u — v — w 
are 1 on one of them and on the rest. Using these facts and Eq |5]) we get 



X= Pr(u, v,w) {H(l — u + w — m + v — z — 0.5) 

U—Uq V — Vq W — Wq 

xH(m - I + 0.5) + H(m -v + z - I + u- w + 0.5)H(l -m - 0.5)} 



m % — I 
K + IT 
K-i + 1 



>i. 



K 



I 

K 



K — i — m 
K 



5l,m+l + 17 Si,- 
K 



(25) 



Using this result in Eq. (24 1 we have 

7 (i) 



1 

%2K 



(.91 + .92) 



(26) 
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where 



9i 



!I2 



K * — ' \ i 

t=0 



i?E(f) EE 

i=0 v 7 Z=o m=C 

i?E(T) EE 



i\ I K — i 



1=0 m=0 



i\ (K — i 



From Eqs. (13) and (27) it follows that 



.91 = 2/i 



1 f2K 



<fi_i, m (iT-* + l), (27) 

(28) 

(29) 



2 V^, 

Thus, we only need to calculate gi- By using Vandermonde's identity and the 
mean value of the hypergeometric function, Eqs. (16) and (17), we obtain 

K 

1 

92 



i=0 
K 



1=0 



i\ I K — i 



I 



-my** 

i=0 

i K 

= -Y 



i=0 
K 



/ 1=0 

K\ (JC 



i 

i-l 



K ^ \ i 



K\( K 



K-i 



K-i 



1 K (2K\ 1 (2K 



K 2\K ) 2\K 



(30) 



case of Eq. JT2l>): 



Using Eqs. Q, (|29j) and Q in Eq. Q, we get the desired result (the 6 = 0.5 

K f2K\ 

S [p = 0.5, oo = l,o fl = -1,0 = 0.5] = 22k( k )• (31) 



E.3 S (p = 0.5, a G = 1, a R = 1, = -0.5) Since, 
([10]) and ([5c| we get 

S =KI^, 

K i K- 



' (1) = iE(f) EE 



2 2K — v % 

i=0 



1=0 m=0 



K-i 
m 



l(a G = l,a R = 1,( 



-0.5, from Eqs. 

(32) 

: -0.5,6 = 0.5). 

(33) 



To calculate I(ac = 1, clr = 1, = —0.5, 6 = 0.5) we look back at the derivation 
of T{ao = 1, clr = 1, ^ = 0, b = 0.5) in the last section (Sec. [5| and notice that 
both cases are similar. The difference lies in that, in this case, the threshold 
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6 = —0.5 < 0. This threshold then makes the nodes in which the sum of the 
updating function equals 0, (namely for which a n j<7 rlj — 0), become active, 
as the sum is larger than the threshold. Therefore, the only cases in which 
damaging one of the inputs of a node changes the state of the target node is 
if either the sum before the damage was and after damage it is below 0, 
or vice versa. Again, since kd — 1, the only possible values for u, v, w and 
z = kd — u — v — w are 1 on one of them and on the rest. Using all this in 
Eq. ([5]) we obtain 



X= Pr(u, v, w) {H(l — u + w — m + v — z + 0.5) 

U—Uq V — Vo W — Wq 

xH(m - I - 0.5) + H(m -v + z-l + u-w- 0.5) HQ, - m + 0.5)} 



m i — I 



K K 
i + 1 



Sl + l.m + 



I K — i — m 
K + 



K 



-fT- 3l+l,m + &l,m — T7 $l,m- 



01.. 



(34) 



Using this result in Eq. (24 1 we get 



= ^K(hi + h 2 -h 3 ), 



(35) 



where 



i=0 v 7 1=0 m=0 v 7 v 7 
i=0 v 7 (=0 m=Q v 7 v 7 
i=0 v 7 i=o m=0 v 7 v 7 



(36) 
(37) 
(38) 



By comparing Eq. ( 14 ) with ( 36 ) and Eq. ( 28 ) with ( 38 ) , we have 



hi = 2/ a 



1 (2K 
2\K 



1 [2K 

■ 92= 2 U 



(39) 
(40) 
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The only term that has not been calculated yet is hi . To do this we use Vander- 



mondc's identity, Eq. (16), with m = n = K 

, I J V I 



i=0 
K 



1=0 



-£(W0(V 
-t(T fK 

i=0 v 

-a 



(41) 



Finally we get the 9 = -0.5 case of Eq. (12) using Eqs. (39), (40), (41) and (351 



in Eq. (32) 



So (p = 0.5, a G = = 1,0 = -0.5) 



K (IK 



2 2K V K 



(42) 



which is, of course, the same as the case 6 = 0.5 because of the symmetry of the 
weights and the probabilities. 
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